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Thermal runaway instability induced by material softening due to shear heating represents a po- 
tential mechanism for mechanical failure of viscoelastic solids. In this work we present a model 
based on a continuum formulation of a viscoelastic material with Arrhenius dependence of viscosity 
on temperature, and investigate the behavior of the thermal runaway phenomenon by analytical 
and numerical methods. Approximate analytical descriptions of the problem reveal that onset of 
thermal runaway instability is controlled by only two dimensionless combinations of physical param- 
eters. Numerical simulations of the model independently verify these analytical results and allow 
a quantitative examination of the complete time evolutions of the shear stress and the spatial dis- 
tributions of temperature and displacement during runaway instability. Thus we find that thermal 
runaway processes may well develop under nonadiabatic conditions. Moreover, nonadiabaticity of 
the unstable runaway mode leads to continuous and extreme localization of the strain and temper- 
ature profiles in space, demonstrating that the thermal runaway process can cause shear banding. 
Examples of time evolutions of the spatial distribution of the shear displacement between the inte- 
rior of the shear band and the essentially nondeforming material outside are presented. Finally, a 
simple relation between evolution of shear stress, displacement, shear-band width and temperature 
rise during runaway instability is given. 
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I. INTRODUCTION 

Many materials are known to exhibit fracture phenom- 
ena that, in apparent contradiction to the expected fail- 
ure behavior usually associated with plastic yielding or 
brittle cracking, are characterized by shear deformation 
localized along a single or a few bands. Ordinarily, one 
would expect brittle fracture to occur by the opening 
of microscopic cracks along the fault plane, while con- 
ventional plastic yielding associated with deformation 
of crystals is known to be pervasive and involve work- 
hardening which seem to be incompatible with the local- 
ization of the deformation in shear bands indicating some 
form of work-softening. This distinct mode of shear fail- 
ure is observed in a variety of materials, including some 
amorphous solids such as polymers [H, [2, Q and bulk 
metallic glasses [1, 0: rocks under high confining 
pressure [8, 9, 10] and crystalline solids deformed rapidly 
in impact experiments [ill . [T^ . 

The formation of shear bands in bulk metallic glasses 
and rocks under high confinement have attracted much 
attention because, in these materials, the mechanism re- 
sponsible for the weakening of the material in the bands 
often initiate instabilities that lead to catastrophic shear 
failure along one dominant band. Accordingly, shear 
banding represents the primary mode of failure in many 
of these material systems. Experimental studies of frac- 
ture surfaces of rock samples subjected to high confining 
pressures demonstrate that the region near the dominant 
band is essentially devoid of microcracks [1, [11] , pointing 
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towards a mechanism leading to viscous sliding without 
cracking. Similarly, studies of fracture surfaces of metal- 
lic glass samples suggest that catastrophic failure is at- 
tributable to a large decrease of the viscosity of the ma- 
terial in the catastrophic shear band [s', ^13] . As a conse- 
quence, these materials appear to fail in a globally brittle, 
but locally ductile, manner. The shear bands in metallic 
glasses, 10-20 nm thick [isj, seem to be accompanied by 
significant local increases in temperature. Vein patterns 
and solidified drops on the fracture surfaces [1, [3, [IB] in- 
dicative of melting of material during catastrophic failure 
as well as indirect experimental measurements of temper- 
ature rises during shear-band operation 6] lend support 
to this view. Similarly, localized shear failure occurring 
in the deeper parts of the Earth's lithosphere, believed 
to be related to earthquakes located several tens of kilo- 
meters below the Earth's surface, may involve substan- 
tial rises in temperature. Geological field observations 
of melted rock in the form of cm thick pseudotachylyte 
layers along shear faults provide evidence for consider- 
able heat dissipation during such failure [1, [13, [H, [11] ■ 
The catastrophic shear failure process thus seems to be 
characterized by spontaneous release of a substantial pro- 
portion of the stored elastic energy as heat in the region 
of the rapidly forming shear band. 

Apart from the observed differences in the qualitative 
macroscopic features of localized shear failure as com- 
pared to ordinary brittle and plastic failure behaviors, 
the circumstances under which shear failure seem to oc- 
cur indicate that this mode of failure may not be ascribed 
to the conventional mechanisms of opening of microscopic 
cracks or crystallographic slip by dislocation motion. For 
instance, the closure of cracks at high confining pressures, 
non-planar crystal structure of minerals and disorder of 
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mineral grain orientations are all factors that inhibit 
these mechanisms from operating in rocks in the Earth's 
interior. In metallic glasses the high degree of structural 
disorder cause dislocations to experience a large number 
of obstacles, reducing their mobility and inhibiting plas- 
tic flow. Because of the absence of these basic weakening 
mechanisms, mantle rocks and bulk metallic glasses can 
demonstrate exceptionally high strengths. However, al- 
though these materials have strengths approaching the 
theoretical shear strength limit at which atomic bonds 
break the fact that shear-band thicknesses are usu- 
ally much larger than interatomic spacing suggests the 
existence of alternative mechanisms responsible for the 
catastrophic shear failure. In metallic glasses it has been 
proposed that granular structure on the microscopic scale 
may be largel y in volved in determining the finite width of 
a shear band [2l| . Discrete element modeling have shown 
that granular packings subjected to compression may fail 
by the formation of shear bands or faults due to dila- 
tancy, predicting a shear-band width of about ten grain 
diameters [23|. However, experimental investigations 
of deformation in a class of shear yielding polymers [ij 
show that, despite of the fact that the molecular and mi- 
crostructural deformation mechanisms are rather differ- 
ent, the phenomenon of shear banding in these materials 
is strikingly similar to the localization of deformation in 
metallic materials. This suggests that the general large- 
scale features of the shear banding phenomenon might be 
appropriately modeled using a continuum formulation as 
long as the grain size is small. 

Several earlier works have investigated the possibility 
of shear failure induced by softening mechanisms facili- 
tating ductile or plastic-like deformation. Mainly two ex- 
planations have been proposed for the observed localiza- 
tion of shear [3, [23j [21, [l^ . The first, based on various 
micromechanical theories developed by Spaepen, Argon, 
Falk and Langer, and others in order to describe plastic- 
ity in amorphous materials [HI, [H, [13, IH, IM H IsH, Ull , 
suggests that material softening due to structural changes 
is a mechanism for strain localization. In this case it is 
usually assumed that the local heat generation during 
deformation is only a secondary effect and is important 
to the evolution of the material in the bands only in the 
later stages of slip. In a recent work, however. Manning, 
Langer and Carlson [s^ proposed that the heat gener- 
ated by plastic deformation is dissipated in the system's 
configurational degrees of freedom and raises an effective 
temperature rather than the usual kinetic temperature. 
Thus it was shown that the effective temperature could 
provide a mechanism for strain localization. 

The second explanation, first proposed by Griggs and 
Baker [34] and later developed by Ogawa f35| in order 
to explain the occurrence of deep- focus earthquakes, in- 
troduces the concept of thermal softening according to 
which the material is weakened primarily due to the effect 
of local heating. Local heating increases the temperature 
which leads to a corresponding decrease in the strongly 
temperature dependent viscosity. Recent theoretical re- 



sults [36| have demonstrated that, even under nonadia- 
batic conditions, the thermal softening mechanism may 
induce a thermal runaway instability exhibiting progres- 
sive strain localization, thus leading to shear-band forma- 
tion and consequent material failure. These results are 
apparently consistent with experimental results on bulk 
metallic glasses [1,[33 showing that shear-band operation 
cannot be fully adiabatic. 

In the present work we expand on the theoretical in- 
vestigations of thermal runaway instability in solids, al- 
ready presented in abbreviated form in Ref. J36]. Our 
model is based on a simple continuum formulation of 
a viscoelastic medium, i.e., the rheology contains both 
viscous and elastic components. The viscous material 
response is supposed to be induced by thermally acti- 
vated processes, yielding a strongly temperature depen- 
dent viscosity. Thus the model accounts for nonelastic 
mechanical responses shown by real materials even be- 
low the conventional elastic limit, such as the well-known 
phenomena of creep or relaxation. The temperature in 
the system is given by the equation for energy conserva- 
tion. The viscoelastic rheology equation, governing the 
mechanical behavior of the material, is then coupled to 
the energy conservation equation through temperature 
dependent viscosity. We shall make the assumption that 
initiation of localized deformation is triggered by small 
(but macroscopic) local heterogeneities or thermal fluc- 
tuations in the otherwise large-scale homogeneous and 
isotropic material. Hence, an increase in strain rate in a 
weaker zone may cause a local temperature rise due to 
viscous dissipation, which weakens the zone even further. 
As a consequence, the local increase in strain rate and 
temperature may amplify strongly because of the effect 
of shear heating-induced thermal softening of the mate- 
rial. Accordingly, catastrophic shear failure may occur 
as a result of thermal runaway instability. 

The paper is organized as follows. In Sec. a simple 
viscoelastic model is introduced and the basic governing 
equations are formulated. In Sec. IIIII analytical meth- 
ods are used to derive the condition for thermal runaway 
to occur in the adiabatic limit and to estimate the re- 
sulting adiabatic temperature rise. Sec. [IV] presents a 
linear analysis for the purpose of determining the con- 
ditions necessary for thermal runaway to occur for the 
general case by taking into account the effects of thermal 
conduction. Numerical solutions to the exact equations 
are presented in Sec. |Vl allowing us to quantify the later 
stages of the thermal runaway process and, in particu- 
lar, the effects of thermal diffusion will be addressed. We 
proceed in Sec. IVII to derive analytically a relatively sim- 
ple relation between the evolution of stress, displacement 
and temperature rise inside the shear band for an adia- 
batic runaway process. The accuracy of this relation is 
then evaluated by comparing it to the numerical results 
for nonadiabatic processes. Finally, we summarize the 
main conclusions in Sec IVIII 
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II. MODEL 

We consider a model consisting of an infinite viscoelas- 
tic slab having a finite width L in the a;-direction; i.e., the 
geometry is that of a solid bounded by a pair of parallel 
infinite planes (see Fig. [T|). We assume that the slab is 
in a condition of simple shear such that the only nonzero 
component of the displacement field is the y-component, 
which we denote by u. Then the shear stress cTxyi^ '^yx)j 
hereafter denoted by a, is constant throughout the slab 
and hence only a function of the time t (see Eq. ([T]) be- 
low). Our purpose is to examine spontaneous modes of 
internal failure in the slab. Therefore, in order to elim- 
inate any additional effects of far-field deformation that 
could either aid or trigger failure, we impose zero velocity 
at the slab's boundaries while we assume that the initial 
shear stress in the slab is ctq. It is thus assumed that the 
shear stress tr has attained a value ao at t = regard- 
less of the slab's loading history which is not considered 
in the present model. Accordingly, the shear stress in 
the slab decreases with time from it's maximum value 
ctq due to relaxation and viscous deformation in the in- 
terior. Our particular model setup with zero velocity 
boundary conditions amounts to searching for the ulti- 
mate conditions for which the slab will fail, that is, if 
it fails at these conditions it is expected to fail earlier 
at any others. The temperature T in the slab initially 
equals a background temperature Ti,g except in the small 
central region having width h and a slightly elevated tem- 
perature Tq. The small thermal perturbation introduced 
in the central region ensures that initiation of localized 
deformation occurs in the neighborhood of the slab's cen- 
ter. The boundaries {x — ±i/2) arc maintained at the 
temperature T^g. 

Assuming inertial effects are negligible, it follows from 
the translational symmetry in the y- and z-directions of 
the present model that the shear stress must satisfy the 
reduced equation for conservation of momentum 



— - 
dx 



(1) 



showing that a is independent of x. Without loss of gen- 
erality the remaining components of the stress tensor are 
regarded as zero. The viscoelastic rheology is represented 
by the Maxwell model [38] for which the strain rate is 
given by 

dv 1 1 da 

" fi{T, <jf^G~dt' ^ ' 

Here v{x,t) is the velocity in the y-direction, G is the 
constant shear modulus and /i(T, a) is the viscosity. The 
first and the last term on the r.h.s. of the equation repre- 
sent the viscous and elastic components of the strain rate, 
respectively. Assuming that the viscosity of the material 
is governed by thermally activated processes, the func- 
tional dependence of the viscosity on temperature and 
shear stress may be approximated as 
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FIG. 1: (Color online) Initial setup of the viscoelastic slab 
model discussed in the text (cross-section in the a;?/- plane). 
The slab is in a state of stress of simple shear with zero veloc- 
ity {v — 0) boundary conditions. The shear stress a, constant 
throughout the slab, initially {t = 0) equals the value ctq. The 
lines show the initial velocity, ^(a;), and temperature, T{x), 
profiles. The shaded region illustrates a small perturbation in 
temperature T = To of width h at the slab center. Elsewhere, 
the background temperature is T = Ti,g. The geometry of 
the strain rate profile concurs with that of the temperature 
profile. 



Here is a pre-exponential constant, n is a constant 
characterizing the dominant creep mechanism, E is the 
activation energy of creep and R = 8.3 JK^^mole"^ is 
the universal gas constant. Thus the viscosity has Ar- 
rhenius dependence on temperature and it is, in general, 
a non-linear function of the shear stress. 

We simplify the mathematical problem by integration 
of Eq. ([2]) to eliminate the velocity from the system of 
equations. Utilizing the zero velocity boundary condi- 
tions, this gives 



Ae «^^T" + 7;^da:= / —dx^O, 4 
G ot J_L dx 



and we thus obtain the equation which governs the time- 
dependence of a: 



da 

'at 



GA 



(5) 



The temperature is determined by the equation for en- 
ergy conservation 



dT 

'at 



dx^ 



1 



dv 
dx 



1 da 
G~dt 



(6) 



where k is the thermal diffusivity and C denotes the heat 
capacity per volume. The last term in Eq. ^ accounts 
for dissipation in the system and thus includes only the 
viscous part of the strain rate. Upon substituting for 
dv/dx the expression in Eq. ([2]), the energy equation be- 
comes 



fi{T,a)=A 



(3) 



dT 
'dt 



d^T A 
'"'d^^C 



(7) 
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Equations ([5]) and ^ with the specified initial and 
boundary conditions constitute a closed set of equations 
for T{x,t) and (j{t). These two equations provide the 
mathematical basis for calculating all the quantities of 
interest during the deformation process since, from the 
solutions to these equations for T and a, one may calcu- 
late the strain rate directly from Eq. We note that, 
because a{t) is independent of x, it follows from Eq. ^ 
that the geometry of the strain rate profile at any instant 
concurs with that of the temperature profile T{x,t). 

From the considerations above it is clear that we can 
infer important information about the deformation pro- 
cesses by studying the temperature rises in the system. 
Useful analytical results can be obtained by invoking ap- 
proximations appropriate for describing the initial stages 
of evolutions of T and a for which temperature rises are 
comparatively small. If instabilities develop, substantial 
temperature rises are possible, but correspondingly large 
decreases in shear stress act against unlimited growth. 
Therefore, to quantify the later stages of the thermal run- 
away processes, we shall take the maximum temperature 
rise, defined as 



where 



AT — T 



Tn 



(8) 



where Tmax is the maximum temperature with respect 
to both time and position, as the appropriate physical 
quantity to study. We solve the system of equations ^ 
and ^ and estimate ATmax using approximate analyti- 
cal (Sec, mil and IIV P and numerical (Sec.|V| methods. 



III. ADIABATIC CASE 

Important insight into the viscoelastic model is gained 
by first examining the limiting case of adiabatic heat- 
ing for which the first term on the r.h.s. of Eq. ([7|) is 
neglected. Then the temperature is determined by the 
reduced equation 



dT A _^ 



dt C 



(9) 



A. Analytical solution 

As aid towards understanding the slab's deformation 
behavior, we first examine the time evolutions of tem- 
perature and shear stress in the initial stages for which 
temperature changes may be considered comparatively 
small. In this limit it is possible to obtain approximate 
analytical solutions to elucidate the deformation prob- 
lem. 

As a first approximation, we insert the initial condition 
for T in Eq. ([5]) and perform the integration over space 
to obtain 



da 
'di 



(10) 



1 - 



h\ ^(Tq, (To) 



(11) 



is a factor which characterizes the initial perturbation. 
At this point, it is convenient to introduce the quantity 
/iQ = A* (To, CTq) and the dimensionless stress and time 



with 



cr ~ t 

(T = , t = 

fJo Tr 



Mo 



2GA, 



(12) 



(13) 



denoting a characteristic time for stress relaxation in the 
system. The simplified dimensionless form of the equa- 
tion is thus 



'dt 



(14) 



This equation can be integrated directly by separation of 
variables and, for the initial condition a = 1, the corre- 
sponding solutions are 



1 



and 



1. 

-< + l 



n > 1 



(15) 



(16) 



Another simplification follows from Taylor expanding 
1/T to first order about the initial temperature Tq in 
the central perturbed zone. Using this approximation, 
the exponential function in Eq. ([9]) may be written 



EiT-Tg) 

— E _ E 2 

Defining the dimensionless temperature 
E{T-To) 

and using Eqs. ([72]) and ^T7\ . we can rewrite Eq. 
dimensionless form as 



de 

'dt 



(17) 



(18) 



on 



(19) 



which must satisfy the initial condition 9 = in the per- 
turbed zone. Here we have introduced the new quantity 



, GCR 



E 



(20) 



having the same dimension as stress. Now, by substitut- 
ing the solutions and for a in Eq. and once 
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again integrating by separation of variables, we find the 
solutions 



= - In 



n = 1 



(21) 



and 




n — 1 ~ 
—t + l 




(22) 



Both solutions (|2ip and ([221) seen to exhibit two dis- 
tinct modes of evolution, depending on the value of ctq. 
Indeed, the solutions are bounded only if < ac- If this 
condition is violated, the temperature grows unlimited at 
the critical times 



tr_r = - In 



71 = 1 



(23) 



and 



n — 1 



n>l. (24) 



The dramatic change in growth of temperature as cto ex- 
ceeds (Tc indicates that, under adiabatic conditions, dc 
plays the role of a critical stress above which thermal 
runaway occurs. 

A few remarks concerning the validity of Eqs. pT|) -([24 |) 
are appropriate here. The main effects of the approxima- 
tions adopted in deriving these equations is that the heat 
production rate becomes too large for large temperature 
rises and that temperatures may increase without limit. 
In contrast, the exact equations always give finite values 
for temperature rises. Nevertheless, solutions (|2ip and 
((22)) are expected to yield reasonable approximations at 
least into the very early stages of development of adia- 
batic thermal runaway. Hence, we also expect that these 
simplified estimates correctly predict the condition for 
adiabatic thermal runaway to occur. That this indeed 
is the case, will be independently verified by numerical 
methods in section IVl 



B. Maximum temperature rise during adiabatic 
thermal runaway 

We now discuss the later stages of adiabatic thermal 
runaway, for which the analytical solutions obtained in 
the previous section are inapplicable. Then, a simple 
analytical estimate of the maximum temperature rise 
(Eq. ([5])) can be made by considering overall energy bal- 
ance as follows. 



Integration of Eq. ([9]) over space gives 

^— = — (T^ / e "Tdx. 
L dt 6 



(25) 



The expression on the r.h.s. can be obtained from 
Eq. (O, giving 



A 



n+l 



e «r dx 



L da 



Combining Eqs. ([25]) and (j26|) we obtain 

L 



d_ 
dt 



Tdx 



2GC 



= 0, 



which may be integrated over time to give 



AT{x,t)dx 



C 2G 



(26) 



(27) 



(28) 



where AT{x, t) = T{x, t) — T{x, 0). Now, it is reasonable 
to assume that the viscous part of the deformation mainly 
occurs within the initially perturbed zone < h/2. 
Then, essentially no heat will be dissipated outside this 
zone, and because diffusion of heat is assumed negligible 
for the present case, the temperature rises in these outer 
regions become vanishingly small. Inside the perturbed 
region, dissipation of heat is expected to be distributed 
monotonously. To understand the reason, recall that in 
this region the initial temperature, and therefore the ini- 
tial viscosity, are uniformly distributed. Accordingly, as 
long as the process is adiabatic, the temperature rise in 
the same region will be approximately uniform, i.e., inde- 
pendent of X. As a result, the integral in Eq. ([^5)) may be 
replaced by a simple integration of uniform temperature 
rise over the region |xl < h/2. Performing the integra- 
tion, we find that the adiabatic maximum temperature 
rise inside the perturbed region is given by 

' — * m n ff 



2GCh- (2^) 
In writing Eq. ()29|) we have made the assumption that 
the runaway process continues until the stress cT{t) in the 
slab is much smaller than ctq. We recognize the quantity 
(Tq /2G as the elastic energy per unit volume stored in the 
slab in its initial state. Eq. (|^5)) reflects the fact that dur- 
ing thermal runaway, and under the assumptions made, 
all the elastic energy in the system spontaneously dissi- 
pates uniformly as heat in the initially perturbed zone. 
As was mentioned earlier, the geometry of the strain rate 
profile concurs at any instant with that of the tempera- 
ture profile. An immediate consequence of the adiabatic 
assumption is therefore that the width of the shear band 
formed during runaway equals the width of the initially 
perturbed zone in which the uniform temperature rise 
occurs. Although the adiabatic approximation provides 
a simple means to obtain important insight into the pro- 
cess, neglect of thermal diffusion may give misleading 
results. The extent to which the adiabatic assumption is 
valid will be investigated by numerical methods in SecfVl 
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IV. CONDITIONS FOR THERMAL RUNAWAY. 
THE GENERAL CASE 



expansion and utilizing Eq. 
linearized equation for S: 



34]), Eq. ([33]) reduces to a 



In this section we attempt to determine the conditions 
necessary for thermal runaway to occur for the general 
case by taking into account the effects of thermal conduc- 
tion. For that purpose, we investigate the stability of the 
coupled equations ([5]) and ([7|) by a linear analysis. Then, 
as an approximation for the initial stages, we insert the 
initial conditions for a and T in Eq. ([5]) and carry out the 
integration over space to obtain the simplified equation 



da _ 
dt ^ 

which has the solution 

(j{t) = CTo 



'2Tr 



1 - 



2Tr 



(30) 



(31) 



Here is the relaxation time given by eq. (jl3p . Using 
this expression, we expand cr""^^ to first order in time, 
yielding 



a{t) 



n+l 



n+1 



1 



(32) 



The expansion is valid for small times t <^ Tr- Equa- 
tion ([7]), determining the temperature rise in the system, 
may now be approximated by substituting for cr""*"^ the 
result obtained in Eq. ((32|) : 



dT 
'dt 



1 



E(_L_i^ 



(33) 



Our aim is to obtain the conditions for which the pertur- 
bation in the central region becomes unstable. Hence, 
in the following discussion we consider only the tem- 
perature inside the central region < h/2. Assume 
that the temperature evolution is continuous and smooth 
and that there exists a steady temperature Tss for which 
dTss/dt = 0. Then Tss satisfies the reduced equation for 
steady flow 



1 - (n+ 1) 



t 

27V 



(34) 



The stability of the steady state solution may be inves- 
tigated by superposing a perturbation S on Tss- The re- 
sulting non-steady temperature T = Tss + <5 must satisfy 
Eq. ((33|) . describing time-dependent flow. For simplicity, 
we restrict the analysis to arbitrarily small 6. In this 
case, since S/Tss < 1, we have 1/T « l/Tss - S/T^s^ and 
we obtain for the exponential term in Eq. (|33p 



1 



ES 



(35) 



where we keep only the leading term in the last expansion 
because ES/RT^^ <^ 1 for arbitrarily small S. Using this 



dS_ 
di 



(36) 



In our search for spontaneous thermal runaway modes 
we are interested in finding the conditions for which the 
initial perturbation Tq instantly starts to increase. If we 
now anticipate that there is a sharp transition bound- 
ary between the stable modes, where Tq instantly drops, 
and the unstable modes, where Tq instantly grows, one 
may expect that Tq plays the role of a steady state right 
at the transition boundary. Accordingly, to obtain the 
conditions for spontaneous runaway modes, one should 
analyze the case in which the initial perturbation Tq it- 
self becomes a steady state and therefore choose Tss ~ Tq 
in Eq. ([36| . This yields the equation 



85 _ d^S 
dt dx'^ 



alE 



^ , 1 - (n+l)— \ 5. 

CflQRT^ V ^ ^2Tr ' 



(37) 



Then, normalizing time by the relaxation time given by 
Eq. (|13p . and introducing new dimensionless variables 



6 = 



E/R' ^ h' 



we recast Eq. p7p into dimensionless form as 



dS 
'dt 



(38) 



(39) 



Here = h'^/n is the characteristic thermal diffusion 
time for diffusion of heat away from the central re- 
gion, and CTc is the characteristic stress introduced earlier 
in ((20|) . Note that two combinations of physical param- 
eters appear in this simplified equation, namely cro/cc 
and Tr/Td- The first combination entered the analysis 
in Sec. IIII Al in which it was interpreted as the factor 
controlling the stability of the system in the adiabatic 
limit. The second combination, which is the ratio of the 
relaxation time to the thermal diffusion time, will modify 
the stability criterion for processes which are not adia- 
batic. The effect of this dimensionless variable on the 
conditions for thermal runaway will be addressed in the 
analysis below. 

Assuming, for simplicity, that S vanish at the "bound- 
aries" X = ±1/2 and that (5 is a positive, even function, 
the general solution to Eq. ([59]) can be represented as a 
sum of particular solutions as 



fm t^ 



(40) 
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with 



X cos (kmX) ■ 



(41) 



Here km = (2m + l)7r denote the frequencies of the per- 
turbation satisfying the required boundary conditions, 
and Bm denote the amphtudes. It follows from the gen- 
eral formula for Fourier cosine coefficients that Bm must 
decrease with increasing m. Furthermore, we note that 
km increases with increasing m. Accordingly, /o is the 
largest term in the series expansion, and it dominates 
the stability of the perturbation in (|40p . To proceed, 
we must analyze the time evolution of the perturbation, 
which is controlled by the argument of the exponential 
in /q. For the case when the coefhcient in the leading 
term in the argument is negative or zero, /o becomes a 
monotonically decreasing function of t. As a result, 5 
decreases with time, implying a stable situation. In the 
opposite case, when the coefficient in the leading term of 
the argument is positive, /o increases with time until it 
reaches a maximum for which dfo/dt = 0. Solving this 
equation, we obtain the characteristic time 



1 



1-(^| -n 



(42) 



above which the perturbation begins to decrease, and the 
corresponding maximum of /o: 

fo{i,imax) ^ Boe-+'^-^> LUeJ J cos(7rJ) . 

(43) 

In Eqs. ((42|) and ((43|) we used that ko — tt. Thus, for cer- 
tain conditions, the solution to this linearized perturbed 
problem predicts only a limited increase in temperature. 
Due to the linear approximation of the exponential func- 
tions in the linear analysis, however, heat is produced at 
a rate which is only a linear function of S. This approx- 
imation is only valid for infinitesimal perturbations and 
severely underestimates the heat production rate as the 
perturbation grows towards finite amplitudes. Neverthe- 
less, we may assume that if the magnitude of 6 becomes 
substantial before t approaches the characteristic time 
tmaxy the perturbation will develop into thermal run- 
away. To proceed, we therefore investigate the growth 
of the quantity 



/o {x, in 



n+l 



/o(i,0) 



We now anticipate that a thermal runaway will develop 
if the quantity above amplifies beyond the characteristic 
factor e, leading to the equation 



(45) 
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FIG. 2: Plots of the critical values of cto/ctc defining the tran- 
sition boundaries between stable deformation modes and un- 
stable thermal runaway modes according to approximate the- 
oretical predictions. The solid curve shows the results of the 
linear stability analysis taking thermal conduction into ac- 
count, as given by the expression in Eq. (|46p . The critical 
values of ao/ac are seen to increase with increasing Tr/rd- 
The dashed line shows the results of the adiabatic analysis 
made in Sec. Illll which completely neglects the effects of ther- 
mal conduction and thus predicts a threshold stress which is 
independent of Tr/rd- 



for which the solution is 




1+4^2^ 
Td 



(46) 



The expression in determines the critical values of 
the dimensionless variables CTo/fc and T^/Td at which 
the transition between stable modes and unstable ther- 
mal runaway modes occurs. Thus, at the transition be- 
tween the two modes, (Jo/(Jc becomes a function of Tr/Td- 
The curve along which the transition occurs divides the 
fo/cc, Tr/Td plane into a "phase diagram", as illustrated 
in Fig. [2] As the solid curve in Fig. [2] shows, the criti- 
cal values of cto/ctc increase with increasing Tr/Td- How- 
ever, for Tr/Td < 10~2 the transition curve stays almost 
vertical, and the critical oq/oc approximately equals its 
lower bound. The lower bound is obtained in the limit 
Tr/Td — > and corresponds to adiabatic processes. In 
this case Eq. (|46|) reduces to 



(47) 



and we thereby recover the result obtained earlier in 
SeclHAl 

Hence, Uc is the smallest stress required for sponta- 
neous thermal runaway to occur. We note, however, that 
the required stress for onset of thermal runaway begins 
to deviate significantly from Oc only when Tr/Td > 1. 
Therefore CTc provides a good estimate for the critical 
stress above which thermal runaway occurs as long as 
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Tr/Td < 1- In this regime, then, onset of thermal runaway 
is to a good approximation independent of the quantities 
which control kinetic processes since these quantities do 
not appear in the expression for Cc (see Eq. ([101) )• The 
results of this linear stability analysis are consistent with 
earlier numerical estimates of initial stages of runaway 



instability in a two-dimensional setup 
timates of CTc were made in Ref. 



[39|. Explicit es- 
and compared to 
typical failure stresses in metallic glasses and rocks un- 
der high confining pressure, showing that the predicted 
values of tXc are within the correct order of magnitude for 
such systems. 

For the situation when t^It^ S> 1 the rate of relaxation 
or creep in the material is very slow compared to the 
process of thermal diffusion, which means that the heat 
produced rapidly flows away from the initially perturbed 
zone. The stress required for thermal runaway to occur 
in this case is therefore much larger than Oc- 

Finally, it should be mentioned that the linear anal- 
ysis presented in this section, and leading to Eq. (jiGj) . 
was based on the greatly simplified equations (|32p and 
([37)1 . The correctness of the obtained results must thus 
be investigated by a proper analysis of the set of the 
fully coupled equations ^ and ([7]). Fortunately, as will 
be demonstrated in Sec. |Vl the results of the linear anal- 
ysis are found to agree extremely well with the results 
obtained from numerical solutions to the complete equa- 
tions. 



V. NUMERICAL APPROACH 

In order to obtain solutions to the complete equations 
([5|) and ([7]) without simplifying assumptions we perform 
numerical simulations. This enables us to capture any 
nonlinear effects in our thermo-mechanical system and 
to study the complete time evolution of T(x^ t), aify and 
the displacement u{x,i). 



A. Dimensionless equations 

The complexity of the problem can be significantly re- 
duced by introducing dimensionless variables. Inspired 
by the analysis in sections IIIII and IIVI we choose the par- 
ticular normalization 



a 



T = 



T 

EjR 



t = 



t 



X 

h 



(48) 



In most situations of interest Tbg <C E/R {E/R is typi- 
cally of the order of 10000 K), and for simplicity we there- 
fore restrict our numerical analysis to the case T^g <C 1. 

The closed set of equations for T and a can then be 
written on dimensionless form as 



dT _ Tr d^T 
di Td dx^ 



(49) 



and 

dd_ 
'dt 



h/L 



/Vi+(1-A)e 



e^o Tfi2_ (-50) 



Similarly, the dimensionless form of the rheology equa- 
tion becomes 



dv 



da 



(51) 



where v = v / {aoh / no) . Displacements are correspond- 
ingly calculated in units of {aoh/fj,o)Tr = aoh/{2ApG). 

The utility of this non-dimensionalization procedure 
lies in that the original problem containing thirteen di- 
mensional parameters now has been reduced to one con- 
taining only eight dimensionless parameters (as can be 
verified upon inspection of the above equations), substan- 
tially reducing the number of necessary numerical runs. 
Moreover, the dimensionless temperature T may now be 
explicitly expressed as a function of the two combina- 
tions of parameters Tr/Td and uq/uc that were suggested 
by the linear analysis in Sec. lIVI as controlling parameters 
for onset of thermal runaway: 



T = /i ( x,i,fo,fbg, -^,n, — , — 

^ Td CTc 



(52) 



To ensure correct numerical results the coupled equa- 
tions (I49p and ([50)1 are solved using a finite-difference 
method with non-uniform mesh and a tailored variable 
time-step. 



B. Temperature rise 



Based on the analysis conducted in Sec. IV Ai we are 
now in a position to investigate the thermal runaway phe- 
nomenon in a self-consistent manner by numerical calcu- 
lations of the complete time evolution of T and a with 
account of heat conduction. Guided by the results ob- 
tained in Sec. Illli we shall begin by studying the max- 
imum temperature rise AT^ax (Eq- ([H])), normalized by 
the adiabatic temperature rise AT^^^ (Eq. ([^O)) ). as a 
function of the governing variables. 

The temperature always attains the maximum value at 
the centre of the slab. Hence, the maximum temperature 
Tmax with respect to both time and position must satisfy 
the equation 



dT 
'dt 



0, 



(53) 



x=0 



from which, by use of Eq. ([52)) . it is easy to show that 
AT„iax = Tmax — Tq is independent of i and t. Then, by 
rewriting AT^ax terms of dimensionless variables, it 
is immediately clear that the maximum temperature rise 
scaled by the adiabatic temperature rise is a function of 
the remaining six dimensionless variables only, i.e., 



ATmax _ r f rf, rf, h Tr (Tq 
^^max \ 1^ Td (Tc 



(54) 
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simplifying the problem even further. 

To examine the behavior of ATmax/ '^T^nax^ '^^ sys- 
tematically varied all six dimcnsionless parameters and 
computed AT^ax from Eqs. ([35]) and ([50]) for each selec- 
tion of fixed parameter values. We present several sets 
of numerical runs in Fig. [3] in terms of contour plots of 
ATjnax / AT^g^^ versus cto/ctc and r,./Trf for different val- 
ues of the remaining parameters n, L/h, Thg and Tq. As 
the series of contour plots show, ATmax/AT^^^ depends 
strongly on the two variables ao/cFc a-nd t^/t^, but is 
rather insensitive to variations in the four remaining pa- 
rameters. In this sense, AT,nax/ AT^^^, when plotted 
against the two combinations of parameters ao/ac and 
Tr/Td, closely resembles a data collapse. It has thus been 
demonstrated that the maximum temperature rise nor- 
malized by the adiabatic temperature rise is, to a quite 
good approximation, a function of the two dimensionless 
variables Oq/uc and t^/t^ alone, as previously suggested 
by the linear analysis (see Sec. llVj) . 

Having identified the controlling variables for the max- 
imum temperature rise in our system, we can now pro- 
ceed to study a representative contour plot, shown in 
Fig. m in more detail. As can be seen, the plot exhibits 
a low-temperature region, corresponding to stable defor- 
mation processes, and a high-temperature region, cor- 
responding to thermal runaway processes. As was pre- 
dicted by the linear analysis in Sec. llVi these two re- 
gions are sharply distinguished by a critical boundary (or 
"transition curve") dividing the tro/cc, Tr/rd plane into a 
"phase diagram" . The location of the critical boundary 
correlates well with the analytical predictions, see Fig. [2] 
for comparison. This verifies that the conditions for spon- 
taneous thermal runaway to occur have been accurately 
determined. A physical explanation of the stability of 
the deformation processes in the low-temperature region 
is that, there, the effects of thermal diffusion and stress 
relaxation dominate over the positive feedback mecha- 
nism. In the high-temperature region the situation is 
exactly the opposite, leading to instability and conse- 
quently thermal runaway. In the following section we 
turn our attention to the high-temperature region, and 
we shall in particular discuss the effects of thermal diffu- 
sion on thermal runaway processes. 



C. Implications of thermal diffusion: Localization 
of deformation and temperature rise during thermal 
runaway 

The high-temperature region exhibited in Fig. 3] is di- 
vided into two domains, as illustrated by red and brown 
colors. An essential observation may be emphasized 
here: The maximum temperature rises in the red domain 
equal the adiabatic maximum temperature rises AT^^^. 
However, the maximum temperature rises in the brown- 
colored region, i.e., in the neighborhood of the critical 
boundary in the high-temperature region, are seen to be 
much larger than AT^^^. Therefore, these two domains 



manifest thermal runaway processes of distinctly differ- 
ent characters, as will be outlined below. 

Let us first address the runaway processes occurring 
in the red-colored domain. An example of the time evo- 
lution of a thermal runaway process in this domain is 
shown in Fig. [S] As can be seen, the process may be di- 
vided into three stages as follows. During the first stage, 
preceding time marker 2 (i.e., t < 0.17rr), a{t) is seen 
to decrease nearly linearly with time while the tempera- 
ture in the slab center Tx=o{t) and the maximum veloc- 
ity Vmax{t) increase gradually. A corresponding gradual 
increase in the spatial distributions of the temperature 
T and the displacement u is noticed in the outer part 
of the initially perturbed zone having width h. How- 
ever, during the second stage, covered by time markers 
2-5 (i.e., in the neighborhood oi t O.lTr^), the shear 
stress a spontaneously drops to zero and is accompa- 
nied by a corresponding explosive rise in the tempera- 
ture Tx=o up to a maximum, whereas one observes an 
explosive increase in the maximum velocity Vmax up to 
a peak value immediately followed by a rapid decrease 
to a much smaller value again. Moreover, during this 
short period of time, a dramatic increase in the displace- 
ment u{x, t) has occurred in the outer part of the initially 
perturbed region, while a major, essentially uniform rise 
in the temperature T(x, t) is observed in the inner part 
of this region (0 < x < 0.4/i). This stage, then, repre- 
sents an extremely rapid thermal runaway process occur- 
ring during a time interval much smaller than the ther- 
mal diffusion time. The effects of thermal diffusion are 
therefore seen to be negligible. Finally, during the third 
stage, succeeding time marker 5 (i.e., t > O.ITt^), v^ax 
decreases smoothly towards zero while essentially no fur- 
ther changes in the remaining quantities can be seen. It 
is worth emphasizing that the rather constant value of 
the temperature T^^o at this last stage stems from the 
thermal diffusion time being too large for any noticeable 
effects of thermal diffusion to be seen within the small 
time interval exhibited in these plots. 

In summary, we observe that thermal runaway pro- 
cesses in the red-colored domain are characterized by a 
rate of heat production which is much greater than the 
rate of heat conduction. Accordingly, the runaway pro- 
cess is well approximated as adiabatic and, as was ex- 
plained in Sec. IIII Bl the elastic energy in the slab is 
therefore dissipated essentially uniformly throughout the 
initially perturbed zone (|a;| < h/2), thus producing a 
maximum temperature rise AT^ax which equals the adi- 
abatic temperature rise AT^^^. Moreover, during the 
runaway process a dramatic increase in the displacement 
occurs, which represents the formation of an adiabatic 
shear band having width of the same order of magnitude 
as the width h of the initially perturbed zone. 

Next, we proceed to analyze the thermal runaway pro- 
cesses occurring in the brown-colored domain, exhibiting 
much larger temperature rises than AT^^^^ . Figure [5] il- 
lustrates an example of a thermal runaway process in this 
domain. Once again, the time evolution of the process 
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FIG. 3: (Color) Contour plots of ATmax/ ^T^^^ versus ao/ac and Tr/rd- Each contour plot represents a set of numerical runs 
for a certain choice of fixed values for the parameters n, L/h, Thg and Tb (values are specified in white boxes). The collection 
of the six contour plots thus constitutes a numerical example of independent variations of all the six dimensionless parameters 
entering Eq. ([54} . A detailed explanation of a contour plot is given in the caption of Fig. [l] 



may be divided into three stages; stage one precedes time 
marker 2 (i.e., t < O.St^), stage two is covered by markers 
2-5 (i.e., t « O.STr), whereas the last stage succeeds time 
marker 5 (i.e., t > O.Sr^). The initial stage is character- 
ized by a relatively large time interval in which the stress 



a decreases approximately linearly with time. During 
the same stage it is seen that the temperature T^^o ^^nd 
the maximum velocity Vmax increase gradually. There 
is basically no change in the spatial distributions of the 
displacement u{x,t) or in the temperature T{x,t). Dur- 
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FIG. 4: (Color) A representative contour plot of ATmai scaled 
by the adiabatic temperature rise AT^^^ as a function of the 
dimensionless variables (to/o"c and Tr/rd. The dark lines rep- 
resent contour lines. The plot exhibits mainly two sharply dis- 
tinguished regions according to small and large temperature 
rises, thus defining a critical boundary dividing the ao/cTc, 
Tr/Td plane into a "phase diagram". The low-temperature 
region, denoted "stable" (blue and green colors), represents 
stable deformation processes. The high-temperature region 
represents unstable thermal runaway processes and is further 
subdivided into two domains, denoted by "self-localizing ther- 
mal runaway" and "adiabatic thermal runaway" (brown and 
red color, respectively), which represents runaway processes 
of different characters (see text). For computational efficiency 
the very late stages of the self-localizing thermal runaway pro- 
cesses have not been fully resolved. The maximum tempera- 
ture rises presented for these processes are therefore underes- 
timated in this plot. See the captions in Figs. [5] and [6] for an 
explanation of the cross and the dot. 



ing the second stage, however, a thermal runaway occurs 
and thus the stress declines spontaneously towards zero. 
Along with this spontaneous stress drop, one observes an 
explosive rise in the central temperature up to a max- 
imum followed by a very rapid decrease. It should be 
emphasized here that, in contrast to what was observed 
in Fig. O the later decrease in Tx=o highlights the im- 
portance of thermal diffusion for this particular type of 
thermal runaway. The maximum velocity is seen to accel- 
erate extremely fast up to a peak value immediately fol- 
lowed by an equally fast decrease to a much smaller value 
again. The dramatic rise in Vmax as one moves from stage 
one to stage two is accompanied by a corresponding ma- 
jor rise in u{x, t) and T{x, t). However, the displacement 
and temperature profiles in this case show strikingly dif- 
ferent features as compared to the same profiles obtained 
in the former case that was shown in Fig. [5l Indeed, in 
this case, large deformation is seen to occur much closer 
to the slab center and an extremely non-uniform rise in 
temperature occurs around the origin. In other words. 
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FIG. 5: An example of a time evolution of an adiabatic ther- 
mal runaway process (see text) corresponding to the location 
of the dot in Fig. |4l For this particular case To ^ 0.0488, 
fbg f» 0.0476, h/L « 0.033, n = 4, Tr/rd ^ 0.027 and 
(7o/(7c ~ 1.83. Panels (a), (b) and (c) display respectively the 
shear stress a in units of ctq, the temperature in the slab center 
Tx^o in units of E/R and the maximum velocity Vmax in units 
of (To/i//io. The two lower panels show a time sequence of spa- 
tial distributions in the vicinity of the initially perturbed zone 
of (d) the displacement it in units of aoh/2ApG and (e) the 
temperature T in units of E/R (for illustrative purposes we 
show profiles only for positive x). The spatial distributions of 
T and u are plotted for six different times corresponding to 
the time markers (black dots) in panels (a)-(c). In all pan- 
els the time t is given in units of the relaxation time Tr, and 
the position x is given in units of the width h of the initially 
perturbed zone. 



the strain and temperature profiles continuously localize 
during the rapid deformation process. Finally, during the 
third stage, Tx^o decreases gradually due to thermal dif- 
fusion and the deformation process terminates as v^ax 
approaches zero. 

In conclusion, we have seen that thermal runaway pro- 
cesses in the brown-colored domain are characterized by 
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FIG. 6: An example of a time evolution of a self-localizing 
thermal runaway (see text) corresponding to the location 
of the cross in Fig H Here To ^ 0.0488, ftg ~ 0.0476, 
h/L « 0.033, n = 4, Tr/ra ~ 0.006 and ao/ac ~ 3.06. Panels 
(a)-(e) are organized in the same manner as outlined in the 
caption of Fig. [S] Note, however, that here the spatial scale 
in panels (d) and (e) is much smaller than the spatial scale of 
the corresponding panels (d) and (e) in Fig. |S] 
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FIG. 7: Two rigid elastic plates sliding past one another on a 
thin viscous layer. 



continuous localization of the temperature and strain 
profiles during deformation, i.e., these runaway processes 
are spatially "self-localizing" . The elastic energy is thus 
dissipated in a zone much narrower than the width of the 
initial perturbation, resulting in maximum temperature 
rises AT^ax which are much larger than the maximum 
adiabatic temperature rises AT^^^. The self- localization 
of these runaway processes clearly arises from the effects 
of thermal diffusion: by diffusion the temperature profile 
initially acquires a peak in the center where the effect 
of the positive feedback mechanism accordingly is max- 
imized. The runaway therefore develops faster in the 
center than in the regions outside and the deformation 
process finally terminates in a highly localized shear band 
with a characteristic width much smaller than the width 
h of the initially perturbed zone. 

Finally, we emphasize that the self-localizing failure 
modes occur at lower values of the shear stress compared 
to the adiabatic modes. Hence, if the material is sub- 
jected to a shear stress large enough to initiate a thermal 
runaway (i.e., a ~ CTc), the failure process is expected to 
be nonadiabatic and to involve a continuous thinning of 
the developing shear band. 



VI. RELATION BETWEEN EVOLUTION OF 
STRESS, DISPLACEMENT AND 
TEMPERATURE RISE 

Thus far the possibility of catastrophic material failure 
by spontaneous thermal runaway has been demonstrated. 
An interesting question yet to be investigated, however, 
is the relation between stress drop, displacement, shear- 
band width and temperature rise during runaway, as now 
discussed. 

For this purpose, let us again consider the viscoelas- 
tic slab shown in Fig. [T] Assume that the viscosity in 
the central region |a;| < h/2 is independent of x and 
much smaller than in the regions outside. In this case 
the system behaves as if two rigid elastic plates slide past 
one another on a thin viscous layer (recall, however, that 
the deforming slab is not far-field driven, but that the 
outer boundaries x = ±L/2 are clamped as internal de- 
formation occurs). This situation is illustrated in Fig. [7| 
where we have defined the initial distribution of the dis- 
placement field u{x,t) to be u{x,0) = 0. In accordance 
with this definition, the shape of the displacement profile 
is as shown in the figure, and we may define a relative 
displacement Ur{t) = u{h/2,t) — u{—h/2,t) between the 
boundaries of the viscous layer. For simplicity, we shall 
study the case when the width of the viscous layer h 
is much smaller than the system size L, i.e., h/L <^ 1. 
Then, denoting the initial stress by ct^, it is clear that 
the maximum possible relative displacement is approxi- 
mately given by 
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We stress here that (Tq is the initial shear stress in the 
slab, whereas u^""^ is the displacement corresponding to 
the final state for which the slab is completely unloaded. 
As the rigid plates continue to slide past one another 
elastic energy stored in the rigid plates is continuously 
dissipated as heat in the viscous layer. If the sliding pro- 
cess occurs in a time short compared to that for thermal 
diffusion, the energy equation for this layer may be writ- 
ten 



dT 



'dt 



(56) 



where 7 denotes the strain inside the viscous central re- 
gion and the term on the r.h.s. of the equation is the 
work dissipated during irreversible viscous flow. It was 
seen in the examples of thermal runaway events in the 
previous section that the shear stress does not remain 
constant during deformation, but drops spontaneously 
towards zero. Hence, a{t) is a strongly varying function 
of time and as a consequence Eq. ((56)) cannot be inte- 
grated over time directly. This apparent complication 
may be circumvented, however, by instead expressing 
both cr and 7 as functions of the displacement. In corre- 
spondence with the adiabatic assumption made here, the 
temperature rise within the viscous layer will be uniform. 
Since the viscosity then is a function of both a uniform 
temperature and shear stress, it follows that the viscos- 
ity and accordingly the strain within the viscous layer is 
uniform. Then, from inspection of Fig. [71 it is seen that 
the uniform strain can be expressed as a function of the 
relative displacement Ur according to the relation 



Ur{t) 



from which we obtain the strain rate 



(i7 
Hi 



1 dUr{t) 

h dt 



(57) 



(58) 



Similarly, in the limit L » h and because the outer 
rigid plates are purely elastic, it can be gleaned from 
Fig. [7] that the stress approximately is given by 



Ur{t) 



(59) 



Substituting these expressions for d-f/dt and (j{t) in 
Eq. ([56]) and using that Ur{t){dur/dt) — l/2{dul/dt), 



Jjjq. i|OU|i aiiu usmg, uiiau ary^jyuar, 

the energy equation takes the form 



C— = - 
dt dt 



G ,A Ur{t) 



(60) 



Next, it follows from Eq. ^ that [G/2L)ur{t) = i((To- 
(t(<)), and the energy equation may thus be written 



C 



dT 



d_ 

dt 



1 



{ao + ait}) 



Ur{t) 
h 



(61) 



thereby eliminating G/L from the equation. The desired 
relation between dynamic quantities can now be obtained 



by direct integration of the energy equation with respect 
to time. If the deformation process terminates at some 
final time tf, then upon integrating from i = to t/, we 
find 



CAT^ 



((To +C7f)ul 



f 



(62) 



Here a/ = a{tf ), AT^ = T{x,tf ) - T{x,0) and = 
Ur{tf ). Not surprisingly, Eq. (|62)) shows that to correctly 
account for decrease in shear stress during deformation, 
one should use the average of the initial and final shear 
stress in calculating the work done by the rigid plates on 
the viscous layer. 

The relation between dynamic quantities in Eq. (j62p 
was obtained under the assumption of adiabatic condi- 
tions. Yet, it was clarified in the previous section that 
thermal runaway processes are greatly affected by ther- 
mal diffusion and therefore not truly adiabatic. As a 
consequence shear was seen to localize to a region of 
width much smaller than the width h of the initially per- 
turbed region. It is therefore of interest to evaluate the 
accuracy of Eq. (|62p for the more general case, including 
self-localizing runaway processes, by numerical methods. 
For this reason proper definitions of the quantities enter- 
ing Eq. ((5^ . valid also for nonadiabatic processes, are 
needed. 

Figure [H] shows that the first stage of stress relaxation 
does not contribute notably to deformation. We define 
the initial shear stress ai associated with the process of 
shear-band formation to be the stress at the instant ti at 
which the curvature of the temporal stress curve d^a/dt^ 
first becomes negative. Thus cr^ is the shear stress in the 
slab right at the onset of thermal runaway. The final 
stress, (Tf, is defined as the stress corresponding to the 
instant tf at which d'^a/dt^ exhibits a maximum, cr/ is 
thus the shear stress in the slab right at the instant at 
which the thermal runaway process terminates. For il- 
lustration, Ui and (7 f are the stresses corresponding to 
time markers 2 and 5 in Fig.[BI^a), respectively. Next, in 
order to define the shear band properly, we must iden- 
tify the point in the displacement profile where a sharp 
transition from very large strain to small strain occurs 
(for example, in the displacement profile corresponding 
to time marker 6 in Fig.[Sfd), this sharp transition occurs 
at the position x/h « 0.01). For this reason, we must 
analyze the second derivative of the displacement profile 
d^ujdx^ at the instant of time when the runaway process 
terminates (essentially no further deformation occurs for 
later times). We define the shear band to be the region 
\x\ < w/2, where x = w/2 corresponds to the position 
at which d^u/dx^ becomes a minimum (by symmetry 
d^u/dx^ becomes a maximum at a; = —w/2,). The final 
relative displacement across the shear band is then given 
by ul — u{w/2,tf) — u{~'w/2,tf). In our numerical cal- 
culations we replace the quantities ctq and h in Eq. ([62]) 
by (Ti and w, respectively. Lastly, the temperature rise 
AT-'' in Eq. 1^ is now taken to be the maximum tem- 
perature rise ATmax occurring in the center of the slab. 
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With all the quantities properly defined, and since the 
controlling variables for onset of thermal runaway have 
been identified, it is now possible to investigate the ro- 
bustness of Eq. ([S^ by numerical analysis. Numerical 
calculations of the quantities entering Eq. ((62|) versus the 
controlling physical parameters cro/cc and Tr/rd, where 
the controlling parameters are varied over the same range 
as in Fig. HI show that 



Hence, according to our model, the simple relation in (j62p 
is applicable even to nonadiabatic situations leading to 
self-localizing processes. 



VII. DISCUSSION AND CONCLUSIONS 

The mechanism of spontaneous thermal runaway in 
viscoelastic solids has been analyzed within the frame- 
work of a highly simplified one-dimensional continuum 
model. As a first approximation the model is only in- 
tended to contain the most essential physics governing 
the thermal runaway phenomenon and certainly neglects 
many important effects encountered in the more compli- 
cated experimental situations. 

Among the most crucial simplifications made in the 
present model is perhaps our neglect of significant 
changes of material properties that may result from the 
very large temperature rises associated with the runaway 
instability. The model does not include corrections due 
to effects of melting, although melting of the material 
near the shear band may be possible. Also, the influence 
of melting on shear-band formation could potentially in- 
crease the rate of deformation even further and cause 
propagation of elastic waves, which would then require 
consideration of inertial effects. For these reasons the 
calculations made within the present model are strictly 
valid only until these changes take place, and the model 
may therefore be assumed to reproduce only qualitatively 
the correct deformation behavior during the later stages 
of the thermal runaway process. However, the inertia 
becomes significant only after the onset of localization 
and/or thermal runaway and hence does not affect the 
conclusions of our study. This statement was also checked 
by additional numerical studies based on the full system 
of equations which include the inertial terms. 

Our approach has thus been to investigate the simplest 
possible model which still has some expectation of rep- 
resenting a real physical situation reasonably well. Al- 
though the ideal character of the model should not be 
ignored, it allows for a quantitative treatment of the de- 
formation problem that hopefully provides valuable in- 
formation about the behavior of the thermal runaway 
failure mechanism. 

A basis for the theory is the assumption that the tem- 
perature dependence of the material's viscosity can be 



described approximately by an Arrhenius expression, and 
that it generally has a nonlinear dependence on the shear 
stress. This represents a simple model accounting for 
thermally activated transitions in the solid, believed to 
be responsible for nonelastic behavior below the yield 
stress. As a consequence, thermal runaway instability 
may occur due to shear heating-induced thermal soften- 
ing of the material. 

In order to determine the conditions necessary for ther- 
mal runaway to occur a theoretical analysis was carried 
out. Neglecting the effect of thermal diffusion, i.e., as- 
suming adiabatic conditions during deformation, approx- 
imate analytical solutions supposed to be valid for the 
early stages of time evolutions of temperature and shear 
stress were found. For this case it was shown that the 
system becomes unstable against shear banding due to 
spontaneous thermal runaway if the initial shear stress cto 
becomes larger than a critical stress CTc (Eq- (|20p ). To in- 
vestigate the effects of thermal diffusion on the instability 
conditions, we subsequently performed a linear analysis 
using rather rough approximations of the complete sys- 
tem of equations. The stability condition obtained within 
the adiabatic approximation was then modified to include 
an additional controlling combination of physical param- 
eters, namely the ratio of the relaxation time to the 
thermal diffusion time r^. According to this analysis Cc 
still provides a good estimate for the critical stress above 
which spontaneous thermal runaway occurs as long as 
Tr/Td < 1. However, when Tr/rd ^ 1, the rate of relax- 
ation or creep in the material is very small compared to 
the process of thermal diffusion, and the stress required 
for thermal runaway to occur in this case is therefore 
much larger. These results were independently verified 
by finite amplitude numerical analysis as demonstrated 
by a series of contour plots of the maximum temperature 
rise ATjnax (Fig- El)- Hence we conclude that initiation 
of spontaneous thermal runaway is controlled by the two 
combinations of parameters ao/cFc and Tr/rd only. 

Numerical investigation of the maximum temperature 
rises and displacement profiles during thermal runaway 
instabilities revealed two potential types of thermal run- 
away processes having distinctly different characters. 
The first, occurring under near adiabatic conditions, is 
characterized by essentially uniform temperature rise and 
strain inside the perturbed zone and is therefore referred 
to as adiabatic thermal runaway. The second, occurring 
under nonadiabatic conditions, is characterized by con- 
tinuous localization of the temperature and strain profiles 
during deformation and is accordingly referred to as self- 
localizing thermal runaway. However, the self-localizing 
failure modes occur at lower values of the shear stress 
compared to the adiabatic modes. In materials subjected 
to increasing loads the actual failure process is therefore 
expected to be nonadiabatic. Thus, if the shear stress in 
the material exceeds a critical value of the order of ctc, 
the material starts to internally disintegrate by unload- 
ing the elastic energy stored in the bulk of the medium 
through accelerated creep along a continuously narrow- 
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ing band. Since creep is a thermally activated process, 
this rapid increase in creep is achieved by local rise in 
temperature. In turn, the accelerated creep prevents the 
material from local cooling due to thermal conduction. 
This opens for the possibility that some materials become 
unstable against macroscopic perturbations (that local- 
ize extremely while developing) before reaching the the- 
oretical shear strength limit at which the material would 
break locally, i.e., at the lattice scale. In this way the 
material may fail by ductile deformation at scales much 
smaller than the deforming sample size, and notably at 
scales much smaller than the characteristic width of an 
initial thermal perturbation, but which still are orders of 
magnitude larger than interatomic spacing. 

Finally, a quantitative relation between evolution of 
stress, deformation and temperature rise was obtained for 
adiabatic shear banding processes by analytical methods. 
Numerical calculations of self-localizing thermal runaway 
processes within our simple viscoelastic model have been 
carried out, showing that the same relation is valid also 
for this particular case. In order to establish the robust- 
ness of Eq. (j62p . however, it would be instructive to com- 
pare this relation to the results of improved model calcu- 
lations of the later stages of thermal runaway processes 
taking into account important changes in material prop- 
erties as commented upon above. 

Recent studies [i^ [4l| showed that the theory devel- 
oped here and in Ref. [36| is applicable to explain the 
generation of intermediate-depth earthquakes. In con- 
trast to this study in which we consider infinitesimal 
perturbations (Ap w 1 and Tq ~ Tbg), John et al. [ilj 
considers large amplitude finite perturbations of the sys- 
tem caused by water influence on rheological properties 
of rocks subjected to differential stresses. Even though 
the finite amplitude perturbations distort the data col- 
lapse (such as presented in Fig. [S]), the representation 
of results as function of the two combinations of param- 
eters (Jo/ac and r,./Trf was proved to be useful. Using 
laboratory derived properties of diabase, the typical rep- 
resentative of lower crust rocks, John et al. [4l'| show that 
the self-localizing thermal runaway can be considered as 
a potential mechanism for deep earthquakes. 

A thorough, quantitative comparison of the theory pre- 
sented here with experiments on bulk metallic glasses 
and polymers at various loading conditions and tem- 
peratures is outside the scope of the present discussion. 
Nevertheless, a very rough estimate of the ratios ctq/ctc 
and Tr/Td for bulk metallic glasses is possible. For bulk 
metallic glasses, typical values are C=1.6xl0^ J m^^^K^^ 
(Ref. 0], G=34 GPa (Ref. [42]), E= 100-400 kJ mole^i 
(Ref. [U, 'ii'l). Assuming a temperature T^g « 620 K 
and infinitely-small amplitude of the perturbation, i.e., 
Ap w 1, we obtain CTc = 0.8-2 GPa. Our study shows, 
however, that the critical value of the stress ctq needed to 
initiate self-localizing thermal runaway may differ signif- 



icantly from <Jc if the ratio Tr/Td is large. The estimate 
of this ratio is much less certain. Only the thermal dif- 
fusivity K has a well-established value of 3x10^^ m^s^^ 
(Ref. [d]). From Ref. [HI the viscosity at this tem- 
perature is inferred to be approximately 10^^-10^^ Pa-s. 
The characteristic width h of the initial perturbation, 
representing a macroscopic perturbation in the material, 
is highly uncertain but a lower bound may be estimated. 
Experimental observations Q of the width w of the zone 
affected by deformation around shear bands give w « 1 
/im. Because of the self-localizing nature of the defor- 
mation (Fig. [SJ;), we expect that h ^ w, and hence we 
assume h — 10 /im - 10 mm (limited by the typical sam- 
ple size). These very rough estimates give a value for 
the dimensionless ratio Tr/Td ranging from 10~^ to 10^ 
and the critical value of ctq necessary for self-localizing 
thermal runaway to occur may be expected to be close 
to the critical stress CTc- It should be noted that our 
estimate of o-q ~ ctc ~ 0.8-2 GPa represents the upper 
limit of critical stress in the system. Modifications of 
our idealized model setup by e.g. increasing the inten- 
sity of the initial perturbation (decreasing Ap) and ac- 
counting for the effect of dynamically building up the 
shear stress in the system would significantly decrease 
the stress required to initiate instability 4l|. Without 
detailed analysis of particular cases, a quantitative com- 
parison of our theory to experimental results on bulk 
metallic glasses is, at present, somewhat difficult as the 
above estimate of the ratio Tr/Td clearly is insufficiently 
constrained. Explicit studies invoking appropriate con- 
stitutive behavior should be undertaken; for example, it 
might be more realistic using a Vogel-Fulcher-Tamman 
or Cohen-Grest fiS] dependence of viscosity on tempera- 
ture instead of the Arrhenius dependence employed here. 
However, the self-localizing thermal runaway mechanism 
appears to be compatible with current experimental data 
and should, in our opinion, be considered as a potential 
mechanism governing instabilities also in materials such 
as bulk metallic glasses and glassy polymers. 

The study [4l| also demonstrated that the theory 
may be applicable to the system described not only for 
thermal, but also for rheological perturbations such as 
variations of activation energy of creep E and/or pre- 
exponential constant A. The rheological model of the 
system (Eq. may also be extended t_oJnclude low tem- 
perature plasticity (Peierls plasticity 13]). 
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